function diff = KFE_firms(param,aux)
    
    %% Generates transition matrix for firms using binomial approximation
    
    ln_m = linspace(log(param.m_l),log(param.m_u),aux.n_m)';

    % Step size
    Delta=(ln_m(2)-ln_m(1));
    dt=Delta^2./param.sigma^2;

    T_y=sparse(aux.n_m,aux.n_m);

    %% Indices 
    ind_im   = 1+1 + (aux.n_m+1).*(0:aux.n_m-2); 
    ind_ip   = 1-1 + (aux.n_m+1).*(1:aux.n_m-1);
    ind_i    = 1   + (aux.n_m+1).*(0:aux.n_m-1);


    %% Transition probability up

    p_tilde=1/2.*(1+ ( ...
        (param.mu-param.sigma.^2./2+ ...
    (1-param.alpha).*(param.zeta+delta_fun(exp(ln_m),param)-eta_fun(exp(ln_m),param))) ... 
        )./param.sigma^2.*Delta);

    % Infinite mean reversion at upper boundary
    p_tilde(end)=0;

    % Numerical check 
    if min(p_tilde)<0 || max(p_tilde)>1 
        error('Probabilities are outside 0/1; increase aux.n_m') 
    end

    % move up and down
    el_i =   zeros(aux.n_m,1);
    el_im =   ones(aux.n_m-1,1).*p_tilde(1:end-1);
    el_ip   =   ones(aux.n_m-1,1).*(1-p_tilde(2:end));

    % Reflecting barrier at ml
    el_i(1) = (1-p_tilde(1));
    el_ip(1) = (1-p_tilde(2)) ;

    % Reflecting barrier at mu
    el_i(end) = p_tilde(end);
    el_im(end) = p_tilde(end-1);

    % Assign values 
    T_y(ind_im) = el_im;
    T_y(ind_ip) = el_ip;
    T_y(ind_i) = el_i;

    clear el* p_tilde

    diff=(T_y-eye(aux.n_m))./dt;

end
